#!/usr/bin/env julia

"""
Reynolds Number Study Example

This example studies the flow behavior at different Reynolds numbers for the 2D 
lid-driven cavity problem, exploring the transition from steady laminar flow to 
complex recirculating patterns.

Run with: julia reynolds_study.jl

Features:
- Flow behavior analysis at Re = 100, 400, 1000, 3200
- Primary and secondary vortex detection
- Flow pattern evolution with increasing Reynolds number
- Critical Reynolds number identification
- Recirculation zone analysis
- Velocity profile comparison at different Re
- Stability and convergence assessment
- Physical flow physics interpretation

"""

# Add the package to the path
push!(LOAD_PATH, "../src")

using LidDrivenCavity
using Printf
using Statistics
using LinearAlgebra

"""
    analyze_flow_patterns(result::LidDrivenResult) -> NamedTuple

Analyze flow patterns and identify key flow features.
"""
function analyze_flow_patterns(result::LidDrivenResult)
    ψ = result.ψ
    ω = result.ω
    u = result.u  
    v = result.v
    x = result.x
    y = result.y
    
    # Remove NaN values for analysis
    ψ_fluid = ψ[.!isnan.(ψ)]
    ω_fluid = ω[.!isnan.(ω)]
    
    if isempty(ψ_fluid)
        return (analysis_failed = true,)
    end
    
    # Stream function extrema (identify vortex centers)
    ψ_min = minimum(ψ_fluid)
    ψ_max = maximum(ψ_fluid)
    ψ_range = ψ_max - ψ_min
    
    # Vorticity extrema
    ω_min = minimum(ω_fluid)
    ω_max = maximum(ω_fluid)
    ω_range = ω_max - ω_min
    
    # Velocity magnitudes
    u_fluid = u[.!isnan.(u)]
    v_fluid = v[.!isnan.(v)]
    u_max = maximum(abs.(u_fluid))
    v_max = maximum(abs.(v_fluid))
    
    # Find primary vortex center (minimum stream function)
    n = length(x)
    primary_vortex_i, primary_vortex_j = 0, 0
    primary_vortex_strength = 0.0
    
    for i in 1:n, j in 1:n
        if !isnan(ψ[i,j]) && ψ[i,j] == ψ_min
            primary_vortex_i, primary_vortex_j = i, j
            primary_vortex_strength = abs(ψ_min)
            break
        end
    end
    
    # Count recirculation zones (simplified detection)
    # Look for local stream function extrema
    local_extrema = 0
    threshold = 0.1 * ψ_range
    
    for i in 2:n-1, j in 2:n-1
        if !isnan(ψ[i,j])
            neighbors = [ψ[i-1,j], ψ[i+1,j], ψ[i,j-1], ψ[i,j+1]]
            neighbors = filter(!isnan, neighbors)
            
            if !isempty(neighbors)
                if (ψ[i,j] > maximum(neighbors) + threshold) || 
                   (ψ[i,j] < minimum(neighbors) - threshold)
                    local_extrema += 1
                end
            end
        end
    end
    
    # Estimate Reynolds number characteristics
    Re = result.options.Re
    characteristic_time = result.solve_time
    
    return (
        analysis_failed = false,
        Re = Re,
        psi_min = ψ_min,
        psi_max = ψ_max, 
        psi_range = ψ_range,
        omega_min = ω_min,
        omega_max = ω_max,
        omega_range = ω_range,
        u_max = u_max,
        v_max = v_max,
        primary_vortex_strength = primary_vortex_strength,
        primary_vortex_location = (primary_vortex_i, primary_vortex_j),
        estimated_vortices = local_extrema,
        solve_time = characteristic_time,
        converged = result.converged,
        iterations = result.iterations
    )
end

"""
    extract_velocity_profiles(result::LidDrivenResult) -> NamedTuple

Extract characteristic velocity profiles for comparison.
"""
function extract_velocity_profiles(result::LidDrivenResult)
    x, y = result.x, result.y
    u, v = result.u, result.v
    n = length(x)
    
    # Centerline profiles
    center_i = div(n, 2)
    center_j = div(n, 2)
    
    # Vertical centerline (u-velocity profile)  
    u_centerline = Float64[]
    y_centerline = Float64[]
    
    for j in 1:n
        if !isnan(u[center_i, j])
            push!(u_centerline, u[center_i, j])
            push!(y_centerline, y[j])
        end
    end
    
    # Horizontal centerline (v-velocity profile)
    v_centerline = Float64[]
    x_centerline = Float64[]
    
    for i in 1:n
        if !isnan(v[i, center_j])
            push!(v_centerline, v[i, center_j])
            push!(x_centerline, x[i])
        end
    end
    
    # Find maximum values and locations
    u_max_val = isempty(u_centerline) ? 0.0 : maximum(u_centerline)
    v_min_val = isempty(v_centerline) ? 0.0 : minimum(v_centerline)
    
    u_max_location = isempty(u_centerline) ? 0.0 : y_centerline[argmax(u_centerline)]
    v_min_location = isempty(v_centerline) ? 0.0 : x_centerline[argmin(v_centerline)]
    
    return (
        u_centerline = u_centerline,
        y_centerline = y_centerline,
        v_centerline = v_centerline,
        x_centerline = x_centerline,
        u_max = u_max_val,
        v_min = v_min_val,  
        u_max_y = u_max_location,
        v_min_x = v_min_location
    )
end

"""
    reynolds_number_study(reynolds_numbers; grid_size=32)

Study flow behavior across different Reynolds numbers.
"""
function reynolds_number_study(reynolds_numbers; grid_size=32)
    println("\n🌪️  Reynolds Number Study")
    println("="^50)
    println("Grid size: $(grid_size)×$(grid_size)")
    println("Reynolds numbers: $(reynolds_numbers)")
    println()
    
    # Storage for results
    results = []
    flow_analyses = []
    velocity_profiles = []
    solver_performance = []
    
    for Re in reynolds_numbers
        println("🔍 Reynolds Number = $(Re)")
        println("-"^30)
        
        # Adjust solver parameters based on Reynolds number
        max_steps = Re <= 400 ? 2000 : 4000
        tolerance = Re <= 1000 ? 1e-6 : 1e-5
        
        options = LidDrivenOptions(
            n = grid_size,
            Re = Re,
            solver = :julia,
            preconditioner = :diagonal,
            max_steps = max_steps,
            tol = tolerance,
            verbose = false,
            save_history = true
        )
        
        print("  Solving... ")
        flush(stdout)
        
        try
            t_start = time()
            result = solve_lid_driven_2d(options)
            wall_time = time() - t_start
            
            if result.converged
                println("✅ Converged ($(result.iterations) iter, $(@sprintf("%.2f", result.solve_time))s)")
                
                # Analyze flow patterns
                flow_analysis = analyze_flow_patterns(result)
                velocity_profile = extract_velocity_profiles(result)
                
                push!(results, result)
                push!(flow_analyses, flow_analysis)
                push!(velocity_profiles, velocity_profile)
                push!(solver_performance, (
                    Re = Re,
                    converged = true,
                    iterations = result.iterations,
                    solve_time = result.solve_time,
                    wall_time = wall_time,
                    residual = result.residual_norm
                ))
                
                # Display key flow characteristics
                println("  Primary vortex strength: $(@sprintf("%.4f", flow_analysis.primary_vortex_strength))")
                println("  Stream function range: [$(@sprintf("%.4f", flow_analysis.psi_min)), $(@sprintf("%.4f", flow_analysis.psi_max))]")
                println("  Vorticity range: [$(@sprintf("%.2f", flow_analysis.omega_min)), $(@sprintf("%.2f", flow_analysis.omega_max))]")
                println("  Max |u|: $(@sprintf("%.4f", flow_analysis.u_max))")
                println("  Max |v|: $(@sprintf("%.4f", flow_analysis.v_max))")
                
                if velocity_profile.u_max > 0
                    println("  Peak u-velocity: $(@sprintf("%.4f", velocity_profile.u_max)) at y=$(@sprintf("%.3f", velocity_profile.u_max_y))")
                end
                if velocity_profile.v_min < 0
                    println("  Peak v-velocity: $(@sprintf("%.4f", velocity_profile.v_min)) at x=$(@sprintf("%.3f", velocity_profile.v_min_x))")
                end
                
            else
                println("❌ Failed to converge")
                push!(solver_performance, (
                    Re = Re,
                    converged = false,
                    iterations = result.iterations,
                    solve_time = result.solve_time,
                    wall_time = wall_time,
                    residual = result.residual_norm
                ))
            end
            
        catch e
            println("❌ Error: $e")
            push!(solver_performance, (
                Re = Re,
                converged = false,
                iterations = 0,
                solve_time = 0.0,
                wall_time = 0.0,
                residual = NaN
            ))
        end
        
        println()
    end
    
    return (
        results = results,
        flow_analyses = flow_analyses,
        velocity_profiles = velocity_profiles,
        solver_performance = solver_performance,
        reynolds_numbers = reynolds_numbers
    )
end

"""
    compare_flow_characteristics(study_results)

Compare flow characteristics across different Reynolds numbers.
"""
function compare_flow_characteristics(study_results)
    flow_analyses = study_results.flow_analyses
    reynolds_numbers = study_results.reynolds_numbers
    solver_performance = study_results.solver_performance
    
    if isempty(flow_analyses)
        println("❌ No successful solutions to compare")
        return
    end
    
    println("📊 Flow Characteristics Comparison")
    println("="^70)
    println("Re      ψ-range    ω-range     Max|u|   Max|v|   Vortex   Iter   Time(s)")
    println("-"^70)
    
    for i in 1:length(flow_analyses)
        if !flow_analyses[i].analysis_failed
            fa = flow_analyses[i]
            sp = solver_performance[i]
            
            println(@sprintf("%4.0f   %8.4f   %8.2f   %7.4f   %7.4f   %6.4f   %4d   %6.2f",
                    fa.Re,
                    fa.psi_range,
                    fa.omega_range, 
                    fa.u_max,
                    fa.v_max,
                    fa.primary_vortex_strength,
                    fa.iterations,
                    fa.solve_time))
        else
            sp = solver_performance[i]
            println(@sprintf("%4.0f   %8s   %8s   %7s   %7s   %6s   %4d   %6.2f",
                    sp.Re, "N/A", "N/A", "N/A", "N/A", "N/A", sp.iterations, sp.solve_time))
        end
    end
    
    # Physical interpretation
    println("\n🌊 Flow Physics Analysis:")
    println("="^35)
    
    successful_analyses = filter(fa -> !fa.analysis_failed, flow_analyses)
    
    if length(successful_analyses) >= 2
        # Trends with Reynolds number
        Re_values = [fa.Re for fa in successful_analyses]
        vortex_strengths = [fa.primary_vortex_strength for fa in successful_analyses]
        omega_ranges = [fa.omega_range for fa in successful_analyses]
        
        println("Flow evolution trends:")
        
        # Vortex strength trend
        if length(Re_values) >= 2
            strength_ratio = vortex_strengths[end] / vortex_strengths[1]
            omega_ratio = omega_ranges[end] / omega_ranges[1]
            
            println(@sprintf("  Vortex strength scaling: %.2fx increase from Re=%.0f to Re=%.0f", 
                    strength_ratio, Re_values[1], Re_values[end]))
            println(@sprintf("  Vorticity range scaling: %.2fx increase", omega_ratio))
            
            if strength_ratio > 2.0
                println("  → Strong vortex intensification with increasing Re")
            elseif strength_ratio > 1.2
                println("  → Moderate vortex intensification with increasing Re")
            else
                println("  → Vortex strength relatively stable")
            end
        end
        
        # Reynolds number regimes
        for (i, fa) in enumerate(successful_analyses)
            Re = fa.Re
            if Re <= 200
                regime = "Steady laminar flow"
            elseif Re <= 800
                regime = "Complex steady state"
            elseif Re <= 2000
                regime = "Transitional regime"
            else
                regime = "Potentially unsteady"
            end
            
            println(@sprintf("  Re = %.0f: %s", Re, regime))
        end
    end
    
    # Solver performance trends
    successful_performance = filter(sp -> sp.converged, solver_performance)
    
    if length(successful_performance) >= 2
        println("\n⚙️  Solver Performance Trends:")
        
        Re_perf = [sp.Re for sp in successful_performance]
        iter_perf = [sp.iterations for sp in successful_performance]
        time_perf = [sp.solve_time for sp in successful_performance]
        
        # Performance scaling
        iter_ratio = iter_perf[end] / iter_perf[1]
        time_ratio = time_perf[end] / time_perf[1]
        
        println(@sprintf("  Iteration scaling: %.2fx from Re=%.0f to Re=%.0f", 
                iter_ratio, Re_perf[1], Re_perf[end]))
        println(@sprintf("  Time scaling: %.2fx", time_ratio))
        
        if iter_ratio > 3.0
            println("  ⚠️  Significant convergence deterioration at high Re")
            println("     Consider: smaller tolerance, more iterations, or adaptive stepping")
        elseif iter_ratio > 1.5
            println("  ℹ️  Moderate convergence impact at high Re (expected)")
        else
            println("  ✅ Good convergence consistency across Re range")
        end
    end
end

"""
    velocity_profile_comparison(study_results)

Compare velocity profiles across different Reynolds numbers.
"""
function velocity_profile_comparison(study_results)
    velocity_profiles = study_results.velocity_profiles
    reynolds_numbers = study_results.reynolds_numbers
    
    if isempty(velocity_profiles)
        return
    end
    
    println("\n📈 Velocity Profile Comparison")
    println("="^40)
    
    println("Centerline velocity characteristics:")
    println("Re      Peak u    u at y    Peak v    v at x")
    println("-"^45)
    
    for i in 1:length(velocity_profiles)
        vp = velocity_profiles[i]
        Re = reynolds_numbers[i]
        
        if !isempty(vp.u_centerline)
            println(@sprintf("%4.0f   %7.4f   %6.3f   %7.4f   %6.3f",
                    Re, vp.u_max, vp.u_max_y, vp.v_min, vp.v_min_x))
        else
            println(@sprintf("%4.0f   %7s   %6s   %7s   %6s", Re, "N/A", "N/A", "N/A", "N/A"))
        end
    end
    
    # Profile shape analysis
    if length(velocity_profiles) >= 2
        println("\n📐 Profile Shape Evolution:")
        
        for i in 1:length(velocity_profiles)
            vp = velocity_profiles[i]
            Re = reynolds_numbers[i]
            
            if !isempty(vp.u_centerline) && length(vp.u_centerline) >= 3
                # Simple profile characteristics
                n_points = length(vp.u_centerline)
                quarter_idx = div(n_points, 4)
                half_idx = div(n_points, 2) 
                three_quarter_idx = 3 * div(n_points, 4)
                
                if quarter_idx > 0 && three_quarter_idx <= n_points
                    # Profile skewness indicator
                    u_quarter = vp.u_centerline[quarter_idx]
                    u_half = vp.u_centerline[half_idx]  
                    u_three_quarter = vp.u_centerline[three_quarter_idx]
                    
                    if u_half > 0.1  # Avoid division by small numbers
                        asymmetry = abs(u_quarter - u_three_quarter) / u_half
                        
                        if asymmetry > 0.5
                            profile_shape = "highly asymmetric"
                        elseif asymmetry > 0.2
                            profile_shape = "moderately asymmetric"
                        else
                            profile_shape = "nearly symmetric"
                        end
                        
                        println(@sprintf("  Re = %.0f: %s profile (asymmetry = %.2f)", 
                                Re, profile_shape, asymmetry))
                    end
                end
            end
        end
    end
end

"""
    critical_reynolds_analysis(study_results)

Identify critical Reynolds numbers and flow transitions.
"""
function critical_reynolds_analysis(study_results)
    flow_analyses = study_results.flow_analyses
    solver_performance = study_results.solver_performance
    reynolds_numbers = study_results.reynolds_numbers
    
    println("\n🎯 Critical Reynolds Number Analysis")
    println("="^45)
    
    successful_analyses = filter(fa -> !fa.analysis_failed, flow_analyses)
    successful_performance = filter(sp -> sp.converged, solver_performance)
    
    if length(successful_analyses) >= 2
        # Look for significant changes in flow characteristics
        Re_values = [fa.Re for fa in successful_analyses]
        vortex_strengths = [fa.primary_vortex_strength for fa in successful_analyses]
        omega_ranges = [fa.omega_range for fa in successful_analyses]
        
        println("Transition indicators:")
        
        # Check for rapid changes between consecutive Reynolds numbers
        for i in 2:length(successful_analyses)
            Re_curr = Re_values[i]
            Re_prev = Re_values[i-1]
            
            strength_change = (vortex_strengths[i] - vortex_strengths[i-1]) / vortex_strengths[i-1]
            omega_change = (omega_ranges[i] - omega_ranges[i-1]) / omega_ranges[i-1]
            
            if abs(strength_change) > 0.3 || abs(omega_change) > 0.5
                println(@sprintf("  Re %.0f → %.0f: Significant flow change detected", Re_prev, Re_curr))
                println(@sprintf("    Vortex strength change: %+.1f%%", strength_change * 100))
                println(@sprintf("    Vorticity range change: %+.1f%%", omega_change * 100))
                
                if Re_curr <= 300 && abs(strength_change) > 0.3
                    println("    → Possible onset of secondary circulation")
                elseif Re_curr <= 1000 && abs(omega_change) > 0.5
                    println("    → Possible approach to critical Reynolds number")
                end
            end
        end
        
        # Convergence difficulty as indicator of approaching instability
        if length(successful_performance) >= 2
            Re_conv = [sp.Re for sp in successful_performance]
            iter_conv = [sp.iterations for sp in successful_performance]
            
            println("\nConvergence trends (instability indicators):")
            for i in 2:length(successful_performance)
                iter_increase = (iter_conv[i] - iter_conv[i-1]) / iter_conv[i-1]
                
                if iter_increase > 1.0  # More than doubled iterations
                    println(@sprintf("  Re %.0f: Significant convergence difficulty (+%.0f%% iterations)", 
                            Re_conv[i], iter_increase * 100))
                    
                    if Re_conv[i] >= 1000
                        println("    → Possible approach to unsteady regime")
                    end
                end
            end
        end
        
        # Known critical Reynolds numbers from literature
        println("\n📚 Literature Critical Reynolds Numbers:")
        println("  Re ≈ 80-120:   First steady secondary vortex")
        println("  Re ≈ 1000:     Multiple secondary vortices") 
        println("  Re ≈ 8000:     Onset of 3D instabilities (in square cavity)")
        println("  Re > 10000:    Transition to unsteady/turbulent flow")
        
        # Compare with current results
        max_Re_tested = maximum(Re_values)
        if max_Re_tested >= 1000
            println("\n🔬 Current study insights:")
            if any(sp.iterations > 2000 for sp in successful_performance if sp.Re >= 1000)
                println("  • High iteration count at Re ≥ 1000 suggests approaching critical regime")
            end
            
            if any(fa.estimated_vortices > 2 for fa in successful_analyses if fa.Re >= 400)
                println("  • Multiple vortex detection indicates complex flow structure")
            end
        end
    else
        println("Insufficient data for critical Reynolds number analysis")
        println("Recommend testing more Reynolds numbers in range [100, 1000]")
    end
end

function main()
    println("🌪️  LidDrivenCavity.jl - Reynolds Number Study")
    println("="^55)
    
    # Display package information
    lid_driven_info()
    println()
    
    println("📋 Study Overview:")
    println("  Purpose: Analyze flow behavior across Reynolds number range")
    println("  Method: Stream function-vorticity formulation")
    println("  Domain: L-shaped cavity with lid-driven boundary condition")
    println()
    
    # Study configuration
    reynolds_sequence = [100.0, 400.0, 1000.0, 3200.0]
    grid_size = 32  # Balance between accuracy and computational cost
    
    println("🔧 Study Configuration:")
    println("  Reynolds numbers: $(reynolds_sequence)")
    println("  Grid resolution: $(grid_size)×$(grid_size)")
    println("  Expected runtime: ~$(length(reynolds_sequence)*2) minutes")
    println()
    
    # Main study
    try
        study_results = reynolds_number_study(reynolds_sequence; grid_size=grid_size)
        
        # Analysis phases
        if !isempty(study_results.flow_analyses)
            compare_flow_characteristics(study_results)
            velocity_profile_comparison(study_results)
            critical_reynolds_analysis(study_results)
            
            println("\n" * "="^60)
            println("📋 Study Summary")
            println("-"^20)
            
            successful_count = count(sp -> sp.converged, study_results.solver_performance)
            total_count = length(study_results.solver_performance)
            
            println("Successful solutions: $(successful_count)/$(total_count)")
            
            if successful_count >= 2
                Re_range = extrema([sp.Re for sp in study_results.solver_performance if sp.converged])
                println("Reynolds number range covered: $(Re_range[1]) - $(Re_range[2])")
                
                # Key findings
                println("\n🎯 Key Findings:")
                println("• Flow complexity increases significantly with Reynolds number")
                println("• Vortex strength and range scale with Re")
                println("• Solver convergence becomes more challenging at high Re")
                println("• Multiple flow regimes observable in tested range")
                
                if Re_range[2] >= 1000
                    println("• High Re solutions approach transitional flow regime")
                end
            end
            
            # Recommendations
            println("\n💡 Recommendations for Future Studies:")
            println("• For Re > 3200: Consider finer grids (64×64) and tighter tolerances")
            println("• For unsteady analysis: Implement time-accurate integration")
            println("• For higher accuracy: Compare with literature benchmark data") 
            println("• For 3D effects: Consider 3D solver for Re > 2000")
            
        else
            println("❌ No successful solutions obtained")
            println("💡 Try with smaller Reynolds numbers or adjusted solver parameters")
        end
        
    catch e
        println("❌ Reynolds number study failed: $e")
        
        if isa(e, OutOfMemoryError)
            println("💡 Try with smaller grid size or fewer Reynolds numbers")
        elseif isa(e, LoadError) || isa(e, MethodError)
            println("💡 Check package dependencies and installation")
        end
    end
end

# Run the Reynolds number study
if abspath(PROGRAM_FILE) == @__FILE__
    try
        main()
        println("\n🎉 Reynolds number study completed successfully!")
        println("   Results provide insights into flow physics across Re range")
    catch e
        println("\n😞 Reynolds number study failed: $e")
        if isa(e, InterruptException)
            println("💡 Study was interrupted - partial results may be available")
        end
    end
end